% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

beta    = beta_p;% 'beta_p ' beta
nu      = nu_p;% 'nu_p   ' nu
alpha   = alpha_p;% 'alpha_p' alpha
delta   = delta_p;% 'delta_p' delta
lambda  = lam_p;% 'lam_p  ' lambda
chi0    = chi0_p;% 'chi0_p ' chi0
chi1    = chi1_p;% 'chi1_p ' chi1
rho_a   = rhoa_p;% 'rhoa_p ' rho_a
sigma_a = siga_p;% 'siga_p ' sigma_a

psi   = psi_p;
n0    = n0_p;
gamma = gam_p;
sigma_rk = (1+nu_p)/(1+alpha_p*nu_p)*siga_p;

%% Generate random shocks

T = 10000;               % simulation periods

rng(1)                   % same random shock  
shocks = normrnd(0,1,T+1,1);  % random shock

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);

%% Simulation

zthre = 0;

kv(1)  = k_ss;
rbv(1) = rb_ss;
nv(1)  = n_ss;
lv(1)  = kv(1)/nv(1);
av(1)  = a_ss;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);
        
        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));
        
        K  = [1 log(kv(t-1)) log(rbv(t-1)) log(nv(t-1)) log(xiv(t))]';
%        x0 = (log(K)).^expmx;
%        x  = prod(x0,2)';
%        if xiv(t) >= 1
            rhseev(t) = exp(sum(K.*eta_rhsee));
            rbv(t)    = exp(sum(K.*eta_rb));
%        else
%            rhseev(t) = exp(x*eta2_rhsee);
%            rbv(t)    = exp(x*eta2_rbar);
%            runv(t) = 1;
%            zthre = max(t,zthre);
%        end
        
        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        % profit and net worth transition depend on bank run
%        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);
%            if pibv(t) > 0
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
%            else
%                nv(t) = chi1*( pibv(t) + nv(t-1) ) + n0;
%            end
%        else
%            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
%            nv(t) = n_ss*(1 - nd/100);
%        end

        lv(t)  = kv(t)/nv(t);
        
        erk_next = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                 + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));
                
        rkhats_next = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next - erk_next) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;

    end
    % End of T period simulation
    
return

%% Pick up bank run events

window = 3;  % how many periods before and after runs to plot

% threshold net worth to define 'normal' times
%thre_n = n_ss;
thre_n = mean(nv(1000:T)) - 2*std(nv(1000:T));

% pick up all runV = 1

norun_num  = 0;
run_num    = 0;
normal_num = 0;
total_num  = 0;

sum_l_norun  = 0;
sum_l_run    = 0;
sum_l_normal = 0;
sum_p_normal = 0;

rkhat_run  = zeros(1);
rkhats_run = zeros(1);
h_run   = zeros(1);
y_run   = zeros(1);
c_run   = zeros(1);
i_run   = zeros(1);
k_run   = zeros(1);
pib_run = zeros(1);
n_run   = zeros(1);
l_run   = zeros(1);
rb_run  = zeros(1);
pr_run  = zeros(1);
shock_run = zeros(1);
xi_run  = zeros(1);
tfp_run = zeros(1);

for t = 1000:T-window
    
    if xiv(t) < 1
        total_num = total_num + 1;
    end
    
    if xiv(t)<1 && nv(t-1) > thre_n
        
        run_num = run_num + 1;
        
        rkhat_run(run_num,1:1+2*window)  = rkhat(t-window:t+window);
        rkhats_run(run_num,1:1+2*window) = rkhats(t-window:t+window);
        h_run(run_num,1:1+2*window)   = hv(t-window:t+window);
        y_run(run_num,1:1+2*window)   = yv(t-window:t+window);
        c_run(run_num,1:1+2*window)   = cv(t-window:t+window);
        i_run(run_num,1:1+2*window)   = iv(t-window:t+window);
        k_run(run_num,1:1+2*window)   = kv(t-window:t+window);
        pib_run(run_num,1:1+2*window) = pibv(t-window:t+window);
        n_run(run_num,1:1+2*window)   = nv(t-window:t+window);
        l_run(run_num,1:1+2*window)   = lv(t-window:t+window);
        rb_run(run_num,1:1+2*window)  = rbv(t-window:t+window);
        pr_run(run_num,1:1+2*window)  = prv(t-window:t+window);
        shock_run(run_num,1:1+2*window) = shocks(t-window:t+window);
        tfp_run(run_num,1:1+2*window) = av(t-window:t+window);
        xi_run(run_num,1:1+2*window)  = xiv(t-window:t+window);
        
        sum_l_run   = sum_l_run + lv(t);
    else
        norun_num   = norun_num + 1;
        sum_l_norun = sum_l_norun + lv(t);
    end
    
%    if lv(t) < thre_l
    if nv(t) > thre_n
        normal_num   = normal_num + 1;
        sum_l_normal = sum_l_normal + lv(t);
        sum_p_normal = sum_p_normal + prv(t);
    end
    
end


avg_l_run   = sum_l_run   / run_num
avg_l_norun = sum_l_norun / norun_num

avg_l_all   = mean(lv(1000:T))
avg_p_all   = total_num / (T-window-1000)

avg_l_normal = sum_l_normal / normal_num
avg_p_normal = sum_p_normal / normal_num

rkhat_run_avg  = mean(rkhat_run,1);
rkhats_run_avg = mean(rkhats_run,1);
h_run_avg   = mean(h_run,1);
y_run_avg   = mean(y_run,1);
c_run_avg   = mean(c_run,1);
i_run_avg   = mean(i_run,1);
k_run_avg   = mean(k_run,1);
pib_run_avg = mean(pib_run,1);
n_run_avg   = mean(n_run,1);
l_run_avg   = mean(l_run,1);
rb_run_avg  = mean(rb_run,1);
pr_run_avg  = mean(pr_run,1);
shock_run_avg = mean(shock_run,1);
xi_run_avg  = mean(xi_run,1);
tfp_run_avg  = mean(tfp_run,1);

return

%% pick up runV = 1 and no run some periods before

norun_period = 10;
run_num2 = 0;

run_run2 = zeros(1);
Rks_run2 = zeros(1);
Rk_run2  = zeros(1);
h_run2   = zeros(1);
y_run2   = zeros(1);
c_run2   = zeros(1);
i_run2   = zeros(1);
k_run2   = zeros(1);
pib_run2 = zeros(1);
n_run2   = zeros(1);
L_run2   = zeros(1);
ERk_run2 = zeros(1);
Rb_run2  = zeros(1);
Pr_run2  = zeros(1);
Exi_run2 = zeros(1);
xi_run2  = zeros(1);

for t = 1000:T-window
    
    if runV(t) == 1
        
        if sum(runV(t-norun_period:t-1)) == 0
        
        run_num2 = run_num2 + 1;
        run_run2(run_num2,1:1+2*window) = runV(t-window:t+window);
        
        Rks_run2(run_num2,1:1+2*window) = RksV(t-window:t+window);
        Rk_run2(run_num2,1:1+2*window)  = RkV(t-window:t+window);
        h_run2(run_num2,1:1+2*window)   = hV(t-window:t+window);
        y_run2(run_num2,1:1+2*window)   = yV(t-window:t+window);
        c_run2(run_num2,1:1+2*window)   = cV(t-window:t+window);
        i_run2(run_num2,1:1+2*window)   = iV(t-window:t+window);
        k_run2(run_num2,1:1+2*window)   = kV(t-window:t+window);
        pib_run2(run_num2,1:1+2*window) = pibV(t-window:t+window);
        n_run2(run_num2,1:1+2*window)   = nV(t-window:t+window);
        L_run2(run_num2,1:1+2*window)   = LV(t-window:t+window);
        ERk_run2(run_num2,1:1+2*window) = ERkV(t-window:t+window);
        Rb_run2(run_num2,1:1+2*window)  = RbV(t-window:t+window);
        Pr_run2(run_num2,1:1+2*window)  = PrV(t-window:t+window);
        Exi_run2(run_num2,1:1+2*window) = shocks(t-window:t+window);
        xi_run2(run_num2,1:1+2*window)  = xiV(t-window:t+window);
        end
    end
end

run_prob_norunbefore = run_num2 / (T-window-1000)

run_run2_avg = mean(run_run2,1);
Rks_run2_avg = mean(Rks_run2,1);
Rk_run2_avg  = mean(Rk_run2,1);
h_run2_avg   = mean(h_run2,1);
y_run2_avg   = mean(y_run2,1);
c_run2_avg   = mean(c_run2,1);
i_run2_avg   = mean(i_run2,1);
k_run2_avg   = mean(k_run2,1);
pib_run2_avg = mean(pib_run2,1);
n_run2_avg   = mean(n_run2,1);
L_run2_avg   = mean(L_run2,1);
ERk_run2_avg = mean(ERk_run2,1);
Rb_run2_avg  = mean(Rb_run2,1);
Pr_run2_avg  = mean(Pr_run2,1);
Exi_run2_avg = mean(Exi_run2,1);
xi_run2_avg  = mean(xi_run2,1);
